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Abstract. In a recent letter [Chatterjee R, Mohanty P K, and Basu A, 2011 J. Stat. 
Mech. L05001], it was claimed that a four-state predator-prey (4SPP) model exhibits 
critical behavior distinct from that of directed percolation (DP). In this Letter, we 
show that the 4SPP in fact belongs to the DP universality class. 
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Absorbing phase transitions have been studied extensively as prototypical emergent 
phenomena of nonequilibrium systems and several universality classes have been 
identified (for a review, see, e. g., Refs. [U El [3l H]). Among them, the directed 
percolation (DP) universality class has been studied thoroughly both numerically and 
analytically. As a result, a so-called 'DP-conjecture' [HI E] has emerged since early 
1980's, asserting that a model should belong to the DP class, if the model exhibits a 
phase transition from an active, fiuctuating phase into a unique absorbing state, if the 
phase transition should be characterized by a one-component order parameter, if no 
quenched disorder is involved, if dynamic rules of the model are short-ranged processes, 
and if the model has neither novel symmetry nor conservation. The validity of the 
conjecture has been supported by numerical evidences. Although the theory suggests 
that the DP class is robust, an experimental realization of the DP class was successful 
only recently [7j. 

Theoretically, it has been questioned if each condition for the conjecture is 
necessary. Indeed, some of conditions are shown to be necessary. It is by now firmly 
established that conditions of no quenched disorder and of short-range processes are 
necessary for the DP critical behavior (for a review of non-DP critical behavior by 
quenched disorder or by long-range interactions, see, e. g., Ref. pj). On the other hand, 
it turns out that certain models that violate some of the above mentioned conditions 
still belong to the DP class. At first, the condition for a unique absorbing state is 
shown to be not necessary. A prototypical model is the pair contact process (PGP) |8] 
which has infinitely many absorbing states (IMAS) but belongs to the DP class. Second, 
the condition of the absence of symmetry or conservation should be applied with care. 
For instance, a parity conservation sometimes changes the universality class as happens 
in branching annihilating random walks with even number of offspring |9l |10], but 
sometimes it has no role as in a model of diffusing particles with reaction dynamics, 
^ and -> [m [12]. 

More notably, the condition of a one-component order parameter is rather tricky 
because this condition requires coarse-graining of microscopic models. In particular, the 
coarse-graining should result in reggeon field theory [13]. A question in this direction 
arises if we have a systematic prescription to conclude that an absorbing phase transition 
of a certain microscopic model is characterized by a one-component order parameter 
field theory. A naive answer would be the number of species in a microscopic model. 
Here, by species is meant a group of particles whose dynamic rules are same irrespective 
of position and time. However, there are counter-examples against this naive answer. 
At first, certain models such as the Ziff-Gulari-Barshad model [H] also exhibit the 
DP criticality, although there are two different species in the model. Second, to make 
things even more complicated, certain single-species microscopic model may have multi- 
component order parameters in field theory. To explain the meaning of the above 
sentence in a concrete way, let us consider single-species models with "hybrid" rules 
such as uA — )• and mA — )■ (m + 1)^4 (n, m > 1) of diffusing A particles [H], which 
show different critical behavior depending on the values of n and m. If m > n, this 
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model belongs to the DP class [llj. However, if m = n, the universality class to which 
this model belongs differs with n. For the case of m = = 1, it still belongs to the DP 
class. For the case of m = n = 2 which is also known as the pair contact process with 
diffusion (PCPD) (for a review see, e. g., Refs. [151 CB]) ^ind for the case of m = n = 3 
which is also known as the the triplet-contact process with diffusion (TCPD) [ITj, these 
models do not share criticality with DP at least in dimensions higher than 1. Regarding 
the number of order parameter components, it was argued by studying the PCPD with 
biased diffusion [19] that the field theory of the PCPD should have at least two 
different order parameter components, which cannot be reduced to a local field theory 
with a single component in the asymptotic regime. Failure of a single-component field 
theory for the PCPD |20] also suggests the same conclusion arrived at in Ref. |18]. To 
summarize, the condition of single component order parameter cannot be answered by 
simply looking at dynamic rules of a model. 

In this context, the conclusion in the recent Letter by Chatterjee et al. [21] is worth 
deep investigation. In that Letter, it was asserted by numerical studies that a four-state 
predator-prey (4SPP) model does not belong to the DP class although the property of 
the model does not seem to deviate much from the conditions of the DP-conjecture. If 
the claim by Chatterjee et al. pT] were right, a theory should be developed why the 
4SPP model does not belong to the DP class, which would deepen our understanding 
of the DP class. In this Letter, we will study the 4SPP model in detail. 

Unfortunately, however, there is an apparent fallacy in Ref. [21]. The 4SPP 
model without diffusion was claimed to have only two absorbing states one of which is 
dynamically unstable. In fact, however, there are infinitely many absorbing states (see 
below). Because of this mistake, Chatterjee et al. [21] chose a wrong order parameter, 
which eventually led them to locate a wrong critical point. In this Letter, we will show 
that the critical point in Ref. [21] is actually underestimated and at the right critical 
point, the 4SPP model exhibits the DP scaling. 

For completeness, we will begin with explaining the 4SPP model. In the 4SPP 
model, there are two species; prey (called A) and predator (called B). At most one 
individual from each species can reside at any lattice point, but two different species 
can stay together at the same site. Hence, each site can have one of four possible 
states; empty, A alone occupied, B alone occupied, and both A and B occupied. A 
can spontaneously generate an offspring to its nearest neighbor, but B cannot branch 
its offspring without consuming A. However, unlike A species, predators also compete 
with each other for space and if two predators are located in a row, both of them can 
be killed together and removed from the system. In this Letter, we restrict ourselves to 
a one-dimensional lattice system with size L. Periodic boundary conditions are always 
assumed. 

To describe the dynamic rules of the model in detail, we introduce a state variable 
Si at site i {1 <i < L) such that Si = sf -\- 2sf, where sf (sf ) is the number of A (B) 
species at site i. Obviously, Si can take one of four possible values, 0, 1, 2, and 3. We 
also assign a bond variable hi to a bond connecting i-th and (i + l)-th sites such that 
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Figure 1. The connection between tlie bond variable hi and the (local) state at sites 
i (circles in left columns) and i + 1 (circles in right columns). Circles in the top 
(below) rows indicate whether A (B) species is present or not. For example, the state 
corresponding to 6i = 11 is that at site i only B is occupied and at site i + I, both A 
and B are occupied. 



hi = 4sj + Sj+i (recall the periodic boundary conditions). The connection between a 
local configuration and a bond variable is depicted in Figure [TJ 

At first, we assume that both A and B are not diffusive. As in Ref. [2T], branching 
occurs in a biased way. In this case, the bond variables changes by the following 
stochastic rules: 

4 — 5 with rate p, 6 — )■ 7 with rate p, 

9 — 7- 10 with rate r, 10 —t- with rate q, 

11 — !■ 1 with rate q, 12—7-13 with rate p, 

. . . , . . I 15 with rate p, 

13 14 with rate r, 14 -> < , ^' 

I 4 with rate q 

15—7-5 with rate q, (1) 

where numbers on the left (right) hand side of the arrows signifies the bond variables 
before (after) a change with corresponding transition rates, p is the transition rate 
of A branching event, q is the transition rate of pair annihilation of B species by 
competition, and r is the transition rate of B branching event by consuming one A. 
Since a configuration change can occur if a bond variable is an element of the following 
set 

S={4, 6, 9, 10, 11, 12, 13, 14, 15}, (2) 

we will refer to a bond with 6j G 5* as an active bond. 

In simulations, we keep a list of active bonds at each time. We choose a time- 
increment dt such a way that 

0<dt< -. (3) 

max(p + q, r) 

Assume that there are N{t) active bonds at time t. Among N{t) active bonds, we choose 
one at random. If hi = 4, 6, or 12, hi increases by 1 with probability pdt (branching 
of A species). In case this change happens, also increases by 4. If 6j = 9 or 13, hi 
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increases by 1 with probability r dt (branching of B species by consuming one A). In 
case this change happens, foj+i also increases by 4. If 6j = 10, 11, or 15, hi decreases 
by 10 with probability qdt (pair annihilation of two predators). In case this change 
happens, (&i-i) also decreases by 8 (2). If hi = 14, 6, increases by 1 with probability 
p dt (branching of A species) or decreases by 10 with probability q dt (pair annihilation 
of two predators). If one of these two events happens, the bond variables of its neighbors 
also changes accordingly. After an attempt of a configuration change, time increases by 
dt/N(t) and the list of active bonds is to be updated. We terminate the above procedure 
if time exceeds the predetermined maximum observation time or the system does not 
have any active bonds. In the following, we choose the transition rates in such a way 
that dt can be chosen to be 1. In all simulations, the initial condition is that = 3 for 
all i. 

Clearly, a configuration without any active bonds is absorbing in the sense that 
no further configuration change can occur. Hence, unlike the claim in Ref. |21], the 
4SPP model without diffusion has infinitely many absorbing states (any configuration 
with only isolated predators is absorbing). Because Chatterjee et al. [21] ignored the 
absorbing nature of isolated predators, the density of predators is wrongfully chosen as 
an order parameter. Since the 4SPP is an IMAS model and the density of predators 
need not be zero in the absorbing state, the critical point reported in Ref. [21] is not 
the correct one. Actually, the density of active bonds defined as 

P«) ^ ^. (4) 

where N{t) is the number of active bonds (see above) at time t and L is the system size, 
can characterize correctly the onset of phase transition. 

At first, we study how p(t)t^ with 6 ~ 0.1594 (DP exponent [22j) behaves near 
criticality. To compare our results with those in Ref. [21], we fix r = 0.9 and q = 0.02 
and the criticality is found by varying p. We simulated the system with size L = 2^° 
and the number of independent runs, for instance, for p = 0.153 81 is 400. Up to the 
observation time tmax ~ 3 x 10^, all simulated systems did not fall into one of absorbing 
states, which minimally guarantees that the finite size effect is not operative. Note that 
the system size and the observation time in this Letter are much larger than those in 
Ref. [21]. In Figure |2]^a), semi- logarithmic plots of the density behavior for various p's 
are depicted. From t = 10^, one can see that the curve for p = 0.153 81 becomes almost 
constant, and the curve for p = 0.153 80 (0.153 82) veers down (up), which indicates 
that the critical point is Pc = 0.153 81(1) with the number in parentheses to be the error 
of the last digit. 

To be more systematic, we also study the effective exponent defined as [b is a certain 
fixed constant) 

_ \npit)/pit/b) 

-m = • (5) 

In the active phase the effective exponent should veer up because p(t) should saturate to 
a finite number. In the absorbing phase, on the other hand, the effective density should 
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Figure 2. (a) Semi-logarithmic plots of p{t)t^ vs t with DP critical exponent 
S = 0.1594 for p = 0.153 75, 0.153 80, 0.153 81, 0.153 82, and 0.153 85 with q = 0.02 and 
r = 0.9. (b) Plots of effective exponents —S{t) as a function of 1/t for the simulation 
results in (a). Since the effective exponent veers down (up) in the absorbing (active) 
phase, the critical point is estimated as pc = 0.153 81(1). Moreover, the effective 
exponent approaches to the DP critical exponent 0.1594 with linear behavior in 1/t. 



veer down because p{t) decays exponentially (or at least faster than the critical decay). 
And only at criticality, 6{t) will approach to the critical exponent 5 as t — t- oo. Note 
that if the leading corrections to scaling takes the 1/t form, that is, if p(t) ~ t~^{l-\-a/t) 
in the asymptotic regime with a constant a, which is the case in many models belonging 
to the DP class, S{t) at criticality approaches to 5 linearly if 6{t) is drawn as a function 
of 1/t. 

We analyze the effective exponent using the same data in Figure |2]^a) with b = 10. 
The resulting plots of the effective exponents are shown in Figure EJ^b). Again we 
conclude that the critical point is 0.153 81(1) and the effective exponent approaches to 
the DP critical exponent 0.1594 linearly in 1/t. 

To make our point stronger, we also study the finite size scaling. To this end, we 
simulated the system at criticality p = 0.153 81 with varying system sizes from L = 2^ 
to L = 2^^. Figure [3]^a) shows the density decay behavior for different system size. Since 
the order parameter takes the form in the scaling regime 

p(t)=F(t/L^)L-/^/^-, (6) 

where z, u^, and (3 are critical exponents and F is a scaling function, plots of p{t)L^^'^^ vs 
t/L^ for different L should collapse into a single curve in the scaling regime. Figure [3]^b) 
shows a nice scaling collapse, using z = 1.58 and /3/z/_l = 0.252 (DP exponents [22]). 
which again confirms that the 4SPP belongs to the DP class. 

Up to now, we assume that both A and B species are not diffusive and we show 
that the 4SPP model with non-diffusive species does belong to the DP class. Now we 
will study if diffusive motion can affect the universal behavior. In Ref. [21], though no 
explicit numerical results are shown, the 4SPP model with diffusion is also claimed not to 
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Figure 3. (a) Plots of ps{t) vs t at criticality {q = 0.02, r = 0.9, p 0.153 81) for 
L = 2^, 2^, 2^", 2", 2^2, and 2^3. (b) Plots of ps{t)L^^''^ vs t/L"" at criticality with 
DP exponents l3/i'± w 0.252 and z « 1.58. As anticipated from Figure[2ja), the data 
collapse becomes good when the time that the finite size effect appears is larger than 
10^ 



belong to the DP class. From the experience of the PCPD and the TCPD, we know that 
diffusion introduced to the model with IMAS may change the universality class. Quite 
interestingly, the exponents reported in Ref. [21] look very similar to those in Ref. [23] (in 
particular see the reported exponents for r = 0.5 in Table II of Ref. [23]), which studies 
a variant of the PCPD. Hence we investigated if Chatterjee et al. indeed observed the 
PCPD behavior from the 4SPP model with diffusion. Since the one-dimensional PCPD 
is known to have strong corrections to scaling and the true scaling behavior has not 
yet been properly observed, Ref. [21] might have observed the intermediated behavior 
of the PCPD, if the 4SPP model with diffusion is PCPD-like. To check this possibility, 
we study the role of biased diffusion, following the suggestion in Ref. [H]. That is, if 
the 4SPP model with biased diffusion shows the same critical behavior as the driven 
PCPD [18], one can conclude that the 4SPP with unbiased diffusion will share criticality 
with the PCPD. 

Let Da {Db) be a transition rate with which A (B) moves to its right nearest 
neighbor. As reasoned above, movement of both species are assumed to be biased. For 
completeness, we fully write the transition rate including diffusion: 

, 5 with rate p, „ \ 7 with rate p 

4-)'<^., „ 6 ^ ' 



2 with rate Da, \ 3 with rate Da, 

8^2 with rate Db, 9 



10 with rate r, 
3 with rate Db, 
10 — !■ with rate q, 11 — )■ 1 with rate q. 



13 with rate p, 
12 <( 9 with rate Da, 13 
6 with rate Db, 



14 with rate r, 
7 with rate Db 
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In practice, we chose transition rates such that dt = 1 can be used near criticahty. 

If Da 7^ and Db = 0, the change of the universality class is not expected, 
because the model still has IMAS and no other symmetry and/or conservation emerges. 
The order parameter defined in Equation (jlj) is again used for = because the 
classification of active bonds does not change. Indeed, as we will see later, the diffusion 
of A does not affect the universal behavior. 

The structure of the absorbing states is significantly affected by diffusion of B. 
If Db 7^ 0, an isolated predator can diffuse around until it meets another A or B. 
Hence the number of absorbing states shrinks to 3; the completely empty configuration, 
a configuration with a single predator without any prey, and the A fully occupied 
configuration without B. 

At first sight, the diffusion of predators seems to have a similar role of the diffusion in 
the PCPD. That is, an isolated predator neither dies out nor branches another offspring 
until it meets another predator or prey. On the other hand, death and branching 
events are mediated by different local states, that is, death requires two predators but 
branching requires one predator and one prey, although in the PCPD both branching and 
annihilation occur by the same local configuration (that is, by a pair) . Actually, whether 
branching and annihilation is mediated by the same local configuration or not seems 
to be an important factor. Unlike the PCPD and the TCPD, the universal behavior of 
models with hybrid rules with m > n remains intact despite of diffusion [HI EU 125] . 
In this regard, the 4SPP with diffusing predators is expected to belong to the DP class 
irrespective of whether the diffusion is biased or not, because this model is similar to 
the model with hybrid-rules of m > n rather than the PCPD. We will show that the 
4SPP model indeed belongs to the DP class. 

Since the diffusion of A is not expected to affect the critical behavior, we study only 
two cases; one is for D^ = 0.1 and Db = and the other is for D^ = and Db = 0.1. q 
and r are fixed to be 0.1 and 0.9, respectively, and p is chosen to be the tuning parameter 
as before. As in the study of the non-diffusive 4SPP model, we study both the critical 
decay and the finite size scaling. The resulting analyses are summarized in Figure ID 
For the model with Da = 0.1 and Db = which has IMAS, the critical point is found to 
be Pc = 0.598 30(5) by studying p(t)t^ with the DP exponent 6 ~ 0.1594. As expected, 
the critical density decay shows DP behavior. The finite size scaling collapse with the 
DP exponents again confirms that the 4SPP with Da = 0.1 and Db = does belong to 
the DP class [see Figure IDj^a)]. 

For the case of Da = and Db = 0.1, a bond is active if a bond variable is an 




15 — )■ 5 with rate q. 



(7) 



Now time-increment dt is chosen to be in the range 



max {p-\- Da + Db, p + q + Da, r + Db) ' 



(8) 
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Figure 4. Numerical Study of critical phenomena for (a) Da — 0.1 and Db = (two 
figures in left column) and for (b) Da — and Db =0.1 (two figures in right column). 
q = 0.1 and r = 0.9 is used. 6 = 0.1594, z = 1.58, and = 0.252 (DP exponents) 

are used. From the two figures in above row, we conclude that Pc = 0.592 30(5) for 
(a) and Pc — 0.431 65(5) for (b). Two figures in lower row shows the finite size scaling 
collapse at corresponding criticality, which concludes that both models belong to the 
DP class. 



element of the following set 

S' = S\J{8} = {4, 6, 8, 9, 10, 11, 12, 13, 14, 15}, (9) 

and now N{t) is the number of these active bonds and the order parameter p{t) is the 
density of such active bonds. In this case, the predator density also can be an order 
parameter. We again analyze the critical decay of the order parameter and the finite 
size scaling. At first, the density decay at pc = 0.431 65(5) is also well described by 
the DP exponent. The finite size scaling confirms that the 4SPP belongs to the DP 
class [see Figure IH^b)]. Hence, we conclude that irrespective of diffusion, the 4SPP does 
belong to the DP class. 

To summarize, we studied the four-state predator-prey model introduced in Ref. [21] 
with and without diffusion. Unlike the conclusion made in Ref. |21], we found that the 
4SPP does belong to the DP class regardless of whether prey and/or predators are 
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diffusive. 
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